{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#  Chapter 6, demo 2\n",
    "\n",
    "### Bayesian data analysis\n",
    "\n",
    "Posterior predictive checking  \n",
    "Binomial example - Testing sequential dependence example"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import preliz as pz\n",
    "\n",
    "import matplotlib.pyplot as plt\n",
    "pz.style.use('preliz-doc')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "# Testing sequential dependence example (Gelman et al p. 163)\n",
    "y = np.array([1, 1, 0, 0, 0, 0, 0, 1, 1, 1, 1, 1, 0, 0, 0, 0, 0, 0, 0, 0])\n",
    "Ty = np.count_nonzero(np.diff(y))\n",
    "\n",
    "# sufficient statistics\n",
    "n = len(y)\n",
    "s = y.sum()\n",
    "\n",
    "nsamp = 10000\n",
    "t = np.random.beta(s+1, n-s+1, size=nsamp)\n",
    "yr = np.random.rand(n, nsamp) < t\n",
    "Tyr = np.count_nonzero(np.diff(yr, axis=0), axis=0)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAABIkAAAH/CAYAAADT8Y0jAAAAOnRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjEwLjYsIGh0dHBzOi8vbWF0cGxvdGxpYi5vcmcvq6yFwwAAAAlwSFlzAAAPYQAAD2EBqD+naQAAXKRJREFUeJzt3Xd4FFX//vF7UyEEQqgKhk5Cb9JEjIA8KKKoFBEpoqKAgIgiIoINlK5IkeaDCA8iIKAgTVEhEEqAhBqKAqGEEkoICYEUdn5/8Nv9ZrObBkk2hPfrurwud2b2zGd2dpfMveecMRmGYQgAAAAAAAD3NRdnFwAAAAAAAADnIyQCAAAAAAAAIREAAAAAAAAIiQAAAAAAACBCIgAAAAAAAIiQCAAAAAAAACIkAgAAAAAAgAiJAAAAAAAAIEIiAAAAAAAAiJAIAHLMjh07FBAQYP1v2LBhzi4pV5w5c8bmuHv06JFtbU+dOtWm7eXLl2db28gazkXe1apVK5tzc79KSkrSd999p+eff1716tWzeU0OHTqUI/vMye8/AAByg5uzCwCAvCqtiytXV1d5eXnpwQcfVK1atdShQwc1atQol6sDAKRn2LBh+u2335xdBnLI+fPn9ffff2vLli06duyYzp8/L8Mw9MADD6hp06Z69dVXVaFCBWeXCQD3HEIiAMiiW7duKTY2VrGxsTp69KiWL1+u3r176/3337fZzt3dXSVKlLA+9vb2zu1SncLV1dXmuH18fJxYDYD70cmTJ+0CooIFC6pQoUKSJDc3/gS+102ePFkrVqywWx4REaGIiAitWLFCs2fPVtOmTZ1QHQDcu/gXEgAyydfXV66urkpOTtbVq1dt1n333Xf6z3/+o3r16lmXNWjQQMHBwblbZB7w4IMP3pfHDSDvCA8Pt3ncsmVLTZ8+Xa6urk6qCDnJ0sM3NjbWuiwhIUHDhw/Xn3/+KZPJ5MTqAODewpxEAJBJP//8s4KDg7Vjxw6tWbPGpreMJP39999OqgwAkNLNmzdtHtesWZOAKB9q0qSJZsyYobCwMO3atUu//vqrihYtal0fGRmp48ePO69AALgH0ZMIAO5A5cqV9fTTT2v+/PnWZdHR0Tbb7NixQz179rQ+fuGFFzR27Fjr42HDhtl0lZ8/f75KlSql6dOna9u2bYqJiVHZsmXVoUMH9e7dO80LnH379mnRokUKDQ1VVFSUbt26pRIlSqhu3brq1KmTHn30UbvnLF++XB9++KH18YABA9S5c2dNnjxZQUFBunHjhgICAvTWW28pMDBQkhQaGqpvv/1We/bs0a1bt1SnTh0NHDhQDRs2tGn7zJkzeuKJJ6yPGzdurAULFlgfX7hwQb/++qsOHDigY8eO6erVq4qJiZG7u7tKlSqlOnXqqEuXLnbtZpfk5GStXr1aa9as0cGDB3X16lUVLFhQVapUUdu2bfXSSy/Jw8PDuv2VK1fUvn17Xbx4UZLk6empX3/9VRUrVrRuM3HiRM2ZM8f6uFevXtbX9/jx41q7dq3279+viIgIXb16VbGxsfL09NSDDz6ohx9+WC+//LKqVatmV6uj94i7u7umT5+uPXv2yM3NTY0bN9bgwYNVqVIlSdIvv/yi+fPn69ixYypYsKAeffRRvffeeypTpoxN247eAz179tS0adP0559/6uLFiypZsqSefPJJ9e/f/46HS+7cuVNLlixRaGioLl26JBcXFz300EMKDAzUq6++ahe25pZWrVopMjLS+vjIkSPavHmz5s6dq/379yspKUn+/v5688039Z///MfmuRm9xyWpR48eCgkJsT7+888/9dBDD6X5/O+++04zZ87UqlWrFBUVpbJly6pr167q0aOHTCaTzp07pylTpigoKEixsbGqUKGCevTooc6dO2fqeLdu3arZs2frwIEDMpvNqlWrlvr06ePw+0GSDMPQ33//rV9++UV79+7VlStX5OHhofLly6t169bq2bOnw/eEo9d17dq1mj9/vo4ePaq4uDib1yIjiYmJWrlypdatW6dDhw4pJiZGBQoU0EMPPaRmzZqpR48eevDBB63bp/7etZg2bZqmTZsmSSpbtqz++uuvTO1fkm7cuKHly5frr7/+0uHDhxUTEyMvLy+VLl1ajRo10osvvujw82thGIZ+/vlnLVq0SMeOHZOHh4cefvhhvfPOO3bPS0xM1NKlS3XgwAEdPnxYly9ftvZcLVq0qKpXr66nn35azz77rFxcbH/rdfRvzmeffabvv/9eK1eu1OnTp1WoUCE9+uijGjx4cJrnYMOGDfrvf/+rw4cPy93dXXXq1FHfvn31wAMPZPi+l6SoqCgtXLhQmzdv1unTp3Xjxg0VK1ZM9evXV/fu3dOcw2/nzp1atGiR9u7dq4sXL8psNqto0aIqWbKkateurfr166t9+/bWfwv79+8vPz8/mzaqVaumJk2aaP369dZlKXsXAQAyRkgEAHcoKSnJ5nHKC5U7sXHjRi1cuFAJCQnWZREREfrqq6905swZjRo1ymZ7wzA0btw4ff/993ZtRUZGKjIyUmvWrFG7du00duxYm+AjtRMnTuj555+3CbrCwsLUp08fTZo0SYZhaOjQoUpOTrau3759u3bv3q158+ZlKdAJCwvTpEmT7JYnJSVZ55JYuXKlBg4cqAEDBmS63cyIiorSgAEDtHfvXrt9h4aGKjQ0VEuXLtWcOXP0wAMPSJKKFSum0aNHq0+fPpJuD2EYOXKkFixYIJPJpIMHD9qcgypVqujdd9+1Pv777781ZcoUu1qSk5P177//6t9//9WyZcv02WefqVOnTunWv27dOv30008ym83WZb///rt27NihJUuWaMGCBfrf//5nXXfz5k399ttv2rVrl3755Rf5+vqm2fb58+f1wgsv2FzgR0ZGau7cudq8ebMWLFiQ7vMdHd/HH3+sZcuW2a07evSojh49qiVLlmjatGlq0qRJptvNKd98842+/fZbm2X79u3TgAEDNGHCBLVv3z7H9n3jxg316NHD5n15/PhxffHFFzp58qRefvllde/eXVeuXLGuP3LkiEaMGKHo6Gi9+eab6ba/cOFCjRo1SoZhWJft2LFDO3bs0CeffKKXX37ZZvu4uDi999572rhxo83yxMREHTx4UAcPHtSSJUs0a9asDO+e9vXXX2vmzJkZvQQOnT17Vm+99ZbdnciSkpJ06NAhHTp0SD/++KO++OILtWvX7o72kZEDBw7o7bfftvlcSFJMTIxiYmJ09OhR+fr6phkSJScna9CgQTahxc2bN/X3338rJCRES5cuVeXKla3r4uLi9Pnnnzts68KFC7pw4YI2btyolStXaubMmXJ3d0+z9ujoaL300ks2Q+8SExP122+/KSQkRL/++quKFStm85xZs2bpq6++slm2efNmBQcHa+DAgWnuy2LDhg0aOnSorl+/blf7unXrtG7dOvXq1UvDhg2zGQK2bNkyffTRRzbvUUm6ePGiLl68qPDwcC1evFht2rSxziuVOiCyHN+BAwesj00mU6YDSQDAbQw3A4AsSkpKUkhIiFavXm1d5uXlddcXkXPnzlVCQoI8PDzsfiFesmSJjh07ZrNs5syZdgGRq6urXRi0evVqjR49Ot19r169WtHR0fLw8LDpsWQ2m/X5559r+PDhSk5OVoECBWyel5SUpAkTJmT6GFNzcXFR4cKF5ePjYzeR7NSpU+3CnLuRmJiovn372rVZqFAhm4uVo0ePql+/fkpMTLQua9Gihbp06WJ9vHPnTi1evFjJycnW10a6PVn5hAkT5Onp6bAGV1dX+fj4qHDhwjbnODk5WZ9++qnOnTuX7jH8+OOPMpvNduchJiZGvXr1sgZEBQsWtFl//vx5zZ07N922f/75Z0VGRsrV1dXuwvOff/5J88I1LWPGjLELiAoUKGDT9rVr1/TWW2/p1KlTWWo7J1gCotSvrSRNmDBBt27dyrF979+/X3v37pXJZLJ77yxcuFCvvfaarly5Ind3d7tzM23aNF27di3d9kePHi3DMOzeF5L05Zdf6siRIzbLhgwZYhcQeXl52Xw3nDt3Tn369LGbny01S0Dk6ekpLy+vdLdNKTExUX369LELiFIfw40bNzR06FDt3LlT0v/dMCB1LycvLy+VKFFCJUqUyHTYeebMGb3++ut2AZGLi4uKFCmSqXluQkNDrQFR6vfW9evXNXXq1HSf7+XlJV9fX7vv9S1btjj8gSCljRs3WgOi1O+rqKgo/fe//7VZtmvXLn399dd27RQoUEBms1nffPNNuvsLCwvTO++8YxMQubi4WEMdi3nz5tl8H5nNZuuPERbu7u5ZvumB2WzWxx9/bHO+WrZs6bTeigBwryIkAoBMeuKJJxQQEKBatWqpR48e1guzsmXLaubMmSpbtuxdtW8ymfTRRx9p9+7d2rp1q+rUqWOzPigoyPr/V65c0axZs2zWDxgwwNobZtSoUTYXdEuWLNE///yT7v7feust7d69W1u2bLE5lujoaCUkJOizzz5TaGio1qxZY3Oxs2fPngwvFFOqXr26Zs6cqeDgYIWHh2vXrl0KCQlRWFiYJk+ebLPt8uXLM91uRlasWKGDBw9aH9epU0fr1q1TaGioduzYYTOkKDw8XL/88ovN84cNG6Zy5cpZH0+YMEFffvmlDh8+bF02YMAA1ahRw+Z5zZo109y5c7Vjxw6Fh4crJCREu3btUmhoqM1wr6SkJK1atSrdY/Dw8LAONVu8eLHNReq5c+fk4+OjH3/8UXv27NGYMWNsnpv6ot+RXr16WWv76KOPbNatXbtWJ06cyLANSTp27Jh+/PFH6+OiRYtq3rx52rNnj/VC0iIuLi7Di8/c4OPjY63xl19+sZnXJCoqyuY854TmzZtr+/btCg0NVatWrazLDcPQ+fPn1aFDB4WEhGjbtm02vXcSEhK0ffv2dNsuVaqUli5dqj179uiPP/6weX5SUpJmz55tfbxlyxab+dXKlSunZcuWKSwsTLt371bXrl2t686dO5dhUOHh4aHx48crNDRUYWFhWrlyZaZCmqVLl+ro0aPWx8WLF9eCBQsUFhambdu2qWXLltZ1ycnJGj9+vKT/u2FA6vfva6+9puDgYAUHBzvs3ebIN998Y/PdVrRoUU2aNElhYWHauXOntm7dqhEjRqh48eLptlOtWjVt2LBBe/bssQvVN2/ebBOOeHl5adKkSfrjjz8UHh6usLAwbd++Xfv27dPy5cttAg9Hd/ZKrVmzZgoODlZYWJhND0fJ9t8USZoxY4ZNLXXq1NHGjRu1Z88ezZkzJ8Mhp+PGjbPpYdu3b1+bHpopX6dp06YpJiZGknTp0iVdvnzZuu7NN99UaGioQkJCtG/fPq1Zs0YjR47Uo48+avcDioXZbNbQoUNtXpMyZcro008/TbdmAIA9QiIAuEsXL17U33//bTME6E488cQT6tmzpzw8POTr66tu3brZrD99+rT1/zdu3KgbN25YH9eqVUsDBw609tR48cUX1bp1a+t6wzBshjukVq5cOb399tvy8PBQsWLF7OYpadKkiV566SW5urqqcuXKdsPLUtaWkfLly6tWrVpatmyZ3nzzTT311FMKDAxUy5Yt7Xo8pe5FcDdS9vySpLFjx1rnFfLx8dHHH39ssz717bO9vLw0btw4a/gWFxenhQsXWtfXr19fb7zxht1+q1evrnLlymnevHnq1auXnnzyST322GNq3bq13TCcjI732WefVevWrWUymVSvXj2VL1/eZv0rr7yihx9+WJL0/PPP24R5Z86cSbftcuXK6YMPPpCXl5c8PDzUs2dPm1tHG4ahTZs2pduGxdq1a20+D++8844eeeQRmUwmubu7q1+/fqpQoYJ1/R9//GHTc8sZ+vfvb62xevXqNnOvSFl7j2eVyWTS559/rqJFi8rNzU1PP/20zfpChQpp5MiR8vLyUuHChe2GVmVU2+DBg62hc7ly5eze60FBQdZwIPXnZOTIkapVq5ak2714RowYYdObJ/XnJLWePXvqueees/YUDAgIsOtZ4si6detsHvfv31+NGzeWyWRSsWLF9OWXX9q8v/ft26ezZ89m2G5mJSYm6o8//rBZNnr0aD3zzDPW/RYrVkw9evSwG66X2tixY+Xn5yeTyaT27dvbDJOKi4uzGeZboEABtWnTRjt27NDgwYP17LPPqkWLFmrevLnefPNNa7Ai3R4mnHqC7pQ8PDw0YcIElShRQq6ururdu7dNT7SU3wk3btywCxtHjRqlBx98UCaTSYGBgXrllVfS3Ne5c+cUFhZmfVy3bl0NHjzY+l6pU6eOzfPj4+OtYaSXl5dN4J0yCPL09FTlypXVvXt3zZ0712FvOEkaP368TchetmxZzZs3T6VLl06zZgCAY8xJBACZ5OvrK1dXV5nNZl29etV6EZyYmKjvv/9enp6eGjx48B23n7L3gCS7X6dThkKpewU1a9bMrr1HHnnEJhhK+at8ak2aNLH5Iz31PBWp54xJr7aMbNu2Tf3797ebs8KRrPRQykjqITWpL8RTSzmvhUWDBg3Uu3dvu15cXl5eGj9+vMPJxVeuXKnhw4fbzWHlSEbHmzK0kW6/JyMiIqyPU54nFxcX+fr6WoewxcfHp9t2kyZN7H6lb9q0qc2FY+ohj2lJ/V779NNP0/1FPyEhQf/8849q1qyZYdujR4/W2rVrHa6bOnWqGjRokKkaU8vK5y+7lStXzqb3XurPX506dWyGamW1tkceecTm8cMPPyx3d3fre/LatWuKiopS6dKl7T4njoLPlM6cOaOrV6/a9LxK6bnnnkv3+WlJ/R5KfQzFihVTQECAzfDRo0eP2k3QfqciIiJsXldvb2+7Ccwzo0yZMqpevbrNsuLFi9sEeyn3c+HCBb3yyiuZ6rVnGIZ1Im9H6tata9PzyNXVVUWLFrVOwp/yO+HkyZM2c86VKFHCbp6lZs2aafr06Q73lbqn3d69ezOcr+rAgQN6/vnn5e3trUaNGlknep85c6a+++47+fn5qXLlyqpevboeffRR1a9f32E7p06d0rx586yPfX19tWDBgrvu3QsA9yt6EgFAJv38888KDg7Wtm3bFBYWptdff91m/YIFC9L9VTcjlomSLVLPPZJyGEDqu7Wkvqh0tCwuLi7Nfae+6Ew9P1DqOR1ShwmpJxtNS2Jiot5///1MBUSSbC5a7lZ6x+/I9evXHfZu6dy5s91cJI0aNbIZimZx6dIlffzxx5kKiKSMjzf1eUj9HsnoPKXH0Xso9bCgzJ63O7mbUOq7A6YlLi5Oly5dcvhfZl9nR7Ly+UvN0bqs1JLR5y/1+qx+/lKfR5PJZBfqWM5tVj8nUvrn7k4v1FPXkZnvuOy8i1XqtlK/PzLLUU+W9N5bo0ePzvSwTin991lm9m1xJ693Snf7mR83bpxNCJScnKwTJ05ow4YNmjp1ql566SW9/PLLNj2pLPbv32/zGrZt25aACADuAj2JAOAOFChQQEOGDNGiRYusv8Zev35dp06dkr+//x21mfrCML1JUQsXLmzzOOVdj9Jalt58Eqn3ndX1mRUWFmb9FVu6PVfK+PHjVbduXXl5eSkxMVG1a9fOln2l5u3tbe2pYzKZMpxHRJLdZMWGYeiTTz6xuyjftGmTfv/9d7Vp08ZmeVBQkE0vgapVq+rLL79UQECAPD09dfz4cbVt2zbTx+Cop1JKd3OeHF3op16W0ZwkFqnfn5ZhVOnJSqCVE1JfPKf3+Uu9ztGF+oULFzK974xem/TuYJUZ0dHRNoGBYRh2vdYs5zb1OS5evHiGEzSnF1JlZmiZIyk/r9Lt77PUwVbq77jU77u7UaRIEZvH58+fv6N2HJ27tF7PxMREm/mg3Nzc9Mknn6hNmzbWY+/SpYv27NmT7ftOfd4d9Wp09O+MRerXvmDBghme+5RDx8qUKaOffvpJhw4d0o4dO3Ts2DFFRETowIED1n9jd+/erWnTptnNN1WvXj2bec2qVKmS7n4BAOkjJAKAbJSTQ1JSqlq1qs3jrVu36r333rNZtm3bNpvHdxpeZaeoqCibx+3atbMZRrJ79+4c23dAQIB27Ngh6fZF7Y8//mg3p09KZrPZLrj43//+p+DgYOtjFxcX67DDjz/+WA0aNLDpzZP6eLt06WIzIXloaOidH1A227FjhwzDsLmITD1HSaVKlTLVlr+/v81QxyFDhqhz585pbu/otU7L2LFjNXbs2Extm1NSz4uSMviUbg+9yc75ce7W9u3bbYZ97d692ybY8vHxUcmSJSXd/pyknOB90qRJdkO9UsrKucsKf39/6/Aj6fb3Wcr335UrV+yGxmXnd1z58uVVsGBB63d6XFycNmzYYDPXW3aLjo62OS8BAQF68cUXrY9jY2MzvAHBnapQoYLNEMSLFy8qIiLCZu6wrVu3pvn81EPLatSoYTN5vSOOwsXq1avbDM+Ljo5Wq1atrEGR5Ts8pbJly9JzCACyEcPNAOAO3Lx5UxMnTrSZ08FkMtlMSJqTWrRoYTMPxYEDBzR16lTdvHlTSUlJWrJkiTZs2GBTW+peLs6Q+tfm4OBg6wX2gQMH7CbUzU6pe+wMGjRIe/bssV6oGIah06dPa8WKFerXr5/dvEPHjx/XxIkTrY8rVqxoc6ei6Ohou1+4Ux/vn3/+qdjYWBmGoW3bttm052wnT57U+PHjdePGDSUmJmr+/Pk2IZHJZFKLFi0y1dZTTz1lExxMmDBB69atsxm+d/nyZf31118aOXKkBgwYkG3HkRt8fX1tel5ERkZqxYoV1vdQyrvW5QVff/21dY6tU6dO6fPPP7dZHxgYaA0HU39OPvroIwUHB9v0qjt//rzWrFmj9957T5999lmO1PzUU0/ZPJ4+fbp27twpwzB05coVDR8+3GZ4b+3atbNtPiLp9qTPqecgGjlypNasWaOEhARJt+dyWrJkiRYtWpQt+/T29rYJaY8fP659+/ZJuh04Dx48ONNDPrOqQIECdhPVjxw5UlFRUTIMQ5s3b9YPP/yQ5vPLlCmjunXrWh/v3r1bX3zxhU1QfvPmTe3bt08zZsxQu3btbG5V/+qrr2rBggX6999/bYbdHj582Pp6S46H5C5fvlwBAQHW/6ZOnZr1FwAAYEVPIgDIpE6dOjmcuNoiMDAw3TkbslOxYsXUt29fm1vGT5s2TTNmzJCrq6vdXDqdO3fOEz2JHn74YXl5eVnDtaNHjyowMFAFCxbU9evX05yANTt07NhRS5YsUXh4uKTbdxLr0qWL3Nzc5O3trevXr9v8ip/yVvbJycn64IMPrBelLi4u+uKLL/Twww/r77//tt7haePGjVqyZIn11/9HH31UJpPJGkRt27ZNTZo0kaenp+Lj43P0eLPKxcVFc+fO1Q8//ODwPdS2bVubXgXpqVKlirp27Wq9+1tMTIwGDRokk8kkHx8fJSQk2PS6a9y4cbYdR26w3O1pzZo11mXDhg3TZ599lmu9CbPi3Llz6tixo03PGAt3d3ebyakDAwPVokULbdy4UdLtAOy1116Tq6urChcurPj4eJv3xgsvvJAjNXfq1Ek//fSTdQLry5cvq3v37g6Pwc3NTR988EG21zBo0CAFBQVZh15duXJFgwcPlouLiwoXLqxr167JMIxsCzkLFSqk+vXrW3sY3rhxQ507d5a3t7d1zqACBQrc1dx36enbt6+2bNli/b4KCQnRY489Zn3NMxp2OGzYMPXs2dP6PTp//nzNnz9fXl5ecnNzswbkjuzfv9/aU8nynZyYmGg34X7KIAoAkDPoSQQAmRQdHa1Lly7pypUrdgFRxYoVNWrUqFytp2/fvurVq5fNslu3btld3Ldr104jR47MxcrSVrhwYb377rs2y8xms65fv24NXnKKh4eHZs2aZXeHnOTkZF29etVuXpmUd5OaOXOm9Rd9Serevbv1VvMjRoywGWI2ZswY652LKlSo4PAcxcfHy9PTM907fuW2rl27qmLFig7fQ1WqVMlyL6/hw4fbDTGzzIWT+iL/Tuetcaa3337bbh4Xy3E1bdo0x+bWuhOWEMNRgDV8+HC7oUKTJk1Sy5YtbZbdunVLV69etXtv5NS58/T01KxZs+zusJX6GAoUKKBx48apUaNG2V7DQw89pP/+9792Q5nMZrNiYmIyPWF/Vnz44Yd24bElIOratavNcNXs1rBhQ4d36Lxx44bc3Nw0ZMgQm+WpQ6MGDRpo8uTJdj0o4+PjrYGahbu7e5pzbVm+k1MHRGXKlNHbb7+dpWMCAGQdPYkA4A54eHjIx8dH/v7+euKJJ9SpUyd5enrmag0mk0kffvih2rVrp0WLFmnXrl2KioqS2WxW8eLFVa9ePXXs2FGPPfZYrtaVkR49eqhkyZL67rvvdPToUXl6eqp27drq06ePmjRpYje3UnYqVaqUFi5cqPXr12vNmjXav3+/rly5IsMw5Ovrq8qVK6thw4Zq3bq19eJ0//79mjFjhrUNPz8/m6DL19dXn376qfVCPD4+XkOHDtXChQvl4uKiYcOGqXz58vrxxx914sQJFSpUSA0bNtTAgQMzPRF0bvD19dXPP/+sadOm6ffff1dUVJRKliypp556Sv37989yrW5ubho9erQ6duyopUuXavfu3YqKilJiYqK8vb1Vvnx51a1bVy1atEh3zpu8qmLFilq0aJG+/vpr7dy5U4mJiapYsaI6dOigbt266dVXX3V2iVYDBw5U9erV9f3331t70tWqVUt9+vRR8+bN7bb39vbWzJkzFRQUpF9//VV79uzRpUuXdOvWLRUpUkQVK1ZU/fr11apVKzVo0CDH6i5TpoyWLl2qX3/9VevWrdOhQ4cUExMjT09P+fn5qVmzZurRo0e2DjNLrVatWvrtt9+0fPly/fnnnzpy5IiuXbumggULqnTp0mrUqJHdsLS7UadOHS1atEjffPONdu3apVu3bqlixYrq2rWrXnzxRfXo0SPb9uVInz59VLlyZX333Xc6fPiw3N3dVb9+fb311lt2Q718fHzsnt+6dWutX79eP/30kzZv3qwTJ04oNjZWnp6eKl26tKpXr65mzZrpP//5j81E5DNmzNC2bdu0a9cunT17VleuXNHNmzdVuHBhVaxYUYGBgerevbvdhOIAgOxnMnLiZxAAAJCnLV++3GbunAEDBmjgwIFOrAhAXjZr1ix99dVX1sd9+/Z12PMIAHBvY7gZAAAAAH399dc28xJZbNmyxW4y/7xwMwQAQPZjuBkAAAAAhYaGaubMmfL19VVAQIBcXV11+vRpnTp1yma7559/XjVr1nRSlQCAnERIBAAAAMAqOjpa27dvd7ju2WefzfUbNQAAcg8hEQAAAAD16dNHFStWVFhYmC5evGiddPrBBx9UvXr11KFDB+udHQEA+RMTVwMAAAAAAICJqwEAAAAAAEBIBAAAAAAAABESAQAAAAAAQIREAAAAAAAAECERAOA+cPToUVWvXl0ff/zxPdFuXnbhwgWNGjVKrVq1Up06ddSmTRvNnj1bZrPZZru4uDg1atRIXbt2dVKl2SuvnOvsfv3PnDmjgIAA9ejR465ry2/nHACA+5GbswsAAKBVq1aKjIy0Webp6alSpUqpUaNGeu2111S1atU7bn/ixIlydXVVnz59bJZ37NhRBw4cyFJbxYoV07Zt29JtN7/asWOHBgwYoGvXrqlMmTIqX768jh49qkmTJunChQsaOXKkdVtvb2/16NFD06dP14YNG9S6dWsnVp4/znVef/3z2jlPadOmTfr+++8VHh6uxMREVaxYUR06dFC3bt3k4nJnv5neSZsbNmxQUFCQDhw4oKioKF29elUFChRQlSpV1LZtW3Xt2lUeHh53c6gAANwVk2EYhrOLAADc3ywhUYUKFVSsWDFJUmxsrCIiIpSUlCQPDw998803atWqVZbb3rVrl7p166YOHTpozJgx1uWGYejVV19VQkKCzfbR0dE6ceKEPDw8VKtWLbv2atasqREjRqTZbn519OhRdenSRWazWWPGjNHTTz8t6faFct++fSVJa9asUcWKFa3PiYmJ0eOPP64yZcpo9erVMplMTqk9P5zrnHr9z5w5oyeeeEKNGzfWggUL7rrOvHLOU5o9e7YmTZokSfLz85OXl5f++ecfmc1mtWrVStOnT89yUHSnbXbt2lWhoaHy8PBQqVKlVLRoUV28eFEXLlyQdPs9N2/ePBUpUuQujxoAgDtkAADgZC1btjT8/f2NZcuW2Sy/ePGi0atXL8Pf399o3LixERcXl+W2Bw0aZPj7+xvbt2/P1PZz5841/P39jW7dumVru/eyW7duGc8//7zh7+9vrFixwm79m2++afj7+xuzZ8+2W/fOO+8Y/v7+xtatW3Oh0qy5V851Tr7+p0+fNvz9/Y3u3btnW7156ZyHhoYaAQEBRrVq1YxVq1ZZlx86dMho1qyZ4e/vb3z33Xe51uayZcuM7du3G4mJiTbLw8LCjMDAQMPf39/49NNPs1QPAADZiTmJAAB5VokSJTR+/Hh5eHjo6tWr2rp1a5aef+XKFW3YsME6bC0zDh8+LEmqVq1atrZ7L1uzZo3Cw8NVu3ZtPf/883br69WrJ0k6duyY3bp27dpJkpYuXZqTJd6Re+Vc32uvf1465zNmzJBhGOrcubOeeeYZ6/Jq1app2LBhkm73CkpKSsqVNjt06KAmTZrI3d3dZnm9evWsz92wYUPmDxAAgGxGSAQAyNNKliypChUqSJIiIiKsywMCAhQQECBJWr9+vbp166aGDRsqICBAZ86ckST98ccfSkpKUmBgYKaHk1iCA0vbjmTU7l9//aWAgAC1b98+zTbi4uLUrFkz1apVy+a48qJ58+ZJkrp37+5wvbe3tyTp6tWrduuaN28uNzc3bdiwQYmJidla18mTJzV//nz17t1bu3btyvLz75Vz7YzX/26OKyfPeVbExcVZg+VOnTrZrX/qqafk7e2tq1evaseOHU5r06JSpUqSpJs3b2bpeQAAZCdCIgBAnmekM33e7Nmz9fbbbysiIsJmTiNJ2rlzpySpTp06mdpPUlKStTdGesFBRu3WrVtX0u2eHWld8M2ePVuXL19Wly5drCFYXnTq1Cnt379fHh4eaU5EbOkx4WjC3QIFCsjf318JCQnat2/fXdVy8+ZNBQUFafTo0WrTpo3atGmjL774Qps3b85STxBLzffCuXbW6383x5Wd5/xuhIeHKykpSZ6enqpRo4bdend3d9WuXVuStHfvXqe1abFnzx5JctguAAC5hbubAQDytIsXL+rkyZOSpPLly9utnzJlikaNGqXOnTvLZDIpOTnZui4sLEzS7clgM+P48eNKSkqSq6ur/P3909wuo3aLFy+usmXLKjIyUocPH7YOB7I4d+6c5s2bp8KFC6t///6Zqs1ZgoKCJEkNGza09lhJ7dq1a5JuH7cjtWvXVnh4uEJDQ9WwYcMs7f/UqVMKCgrSpk2bFBISYg0sXFxcVLduXQUGBiowMNB6YZ5Z98q5dtbrf7fHlZV9zpw5U5s2bcpUXSmNHDky3UDF8r3x4IMPys3N8Z+8fn5+2rZtm3XbjGR3m7du3dLFixf1559/atKkSfLy8tJ7772XqVoAAMgJhEQAgDzr8uXLGjp0qBITE+Xj46NmzZrZbfPSSy/pxRdftD62XLgZhqHz589Luj1kLTMsw4/Kly+vAgUKONwms+3WqVNHkZGROnDggN0F9ldffaWEhAQNGDDApudTXhQaGipJ2rp1a7o9biTHIZ70f69TZGRkhvtLSEhQSEiIgoKCFBQUZDOMqWjRomrdurUCAwP12GOP3dVrd6+c69x+/VO6m+PKyj4jIiKsx5kVsbGx6a6PiYmRJPn4+KS5jeUuYpagLSPZ1ea8efPs7pTXunVrDRo0KN3QEgCAnEZIBADIM2bNmmWd7DY2NlYRERFKSkqSu7u7Ro0a5bAnxXPPPeewrWvXrll7FRUtWjRT+8/MRMaZbbdu3bpau3atDh48aLP8wIEDWrVqlcqUKaNXXnklU3U506FDhyTdfk28vLwcbrNv3z4lJyen+bpZLqijo6Mz3N+IESO0cuVKSZLJZFLNmjX1+OOPKzAwUHXr1s3yrcrTcq+c69x+/VO6m+PKyj7Hjh2rsWPHZqm2zEhISJAku0miU7IM0cvsPEDZ1Wbp0qXVoEEDJScn6+zZs7p06ZJ27Nih1atXq3LlynJ1dc1UPQAAZDdCIgBAnhEREWHtOeLu7q6SJUuqYcOGeu2111S9enWHz6lcubLD5ZaLOUtbmXHkyBFJ6c9Rk9l2LcOfDhw4YLN83LhxMgxD77zzjjw9PTNVl7MYhmGdBHzatGny8/Oz2+bKlStq1qyZ3NzcrPPYpGbpqZOZC/GUQ3gMw1BSUpISEhKUkJCgW7duZVtIdC+ca2e8/indzXHd6T6zk6W29OarskysnVZvspxqs23btmrbtq318d69e/Xxxx9r5syZunr1qj777LNM1QMAQHYjJAIA5BljxoxRhw4dsvSctHpXpOz5ERsbm+7wEAtLcJBe75LMtlurVi25ubnp2LFjunHjhgoWLKgNGzYoJCRENWvWTPeuUdLtgCAxMdGpQVJcXJz1gjet4VZbtmyRYRhq3LixChYs6HAby123fH19M9znF198oS5dumjTpk3atGmTwsPDdfToUf33v/9VoUKF9Mgjj+ixxx5TYGCgypQpc2cHprx1rtPijNc/pbs5rjvdZ3aynC/LEDFHLEPCLEPEnNGmdLvX1uzZs9W6dWstWbJEb775psqWLZvp5wMAkF0IiQAA+ZKHh4e8vb0VFxenmJiYDEOiS5cu6dKlS5LSDw4y267lDk/h4eE6dOiQ6tSpo4kTJ0qSPvjgA5lMJpvthw0bpj179ujTTz/V+PHjdfToUQ0bNkzdu3fXvn37NHnyZOskyk2bNtXw4cNtepZYnv/JJ59o7NixOn78uMqXL6+hQ4cqMDAw4xfMAUtvCZPJlGaviDVr1kiS2rVrl2Y7lgvqzMzJ4+Lionr16qlevXoaNGiQLl++bJ2fKDg4WBs2bNCGDRskSVWqVLFOXN2wYcNM9xhz9rnOLGe8/indzXFlZZ85NXG1ZY6mc+fOKTk52eFE06dPn7bZNiM50aZF6dKlVb16de3du1eHDx8mJAIAOAUhEQAg36pevbp27typY8eOqVy5culua5mjpmjRonrggQeypV3LHZ4OHjyo8PBwnThxQi1btlSTJk0cbn/58mV99NFH6tu3r/z8/FSsWDHt27dP3bp1U9OmTTVhwgQZhqFp06apV69eWrt2rc1tz69cuaLhw4frrbfeUokSJbRgwQK99dZb+uWXX1SlSpV0j8mRIkWKyMXFRWaz2dqTJKUTJ05o06ZNKlmypJ599tk027Hcav5Obu1dvHhxvfDCC3rhhRd069YthYWFafPmzdq0aZMOHz6sf//9V3PnztX333/vcGJzR/LCuc6MvPD63+lxZWWfOTVxdY0aNeTu7q6EhASFh4erTp06NuuTkpK0f/9+SUpzqF5utJmSZQ6sW7duZfm5AABkh+wZ2A8AQB708MMPS7KfU8URS3CQmTsLZbZdy0Xitm3bNG3aNLm6uur9999Pc/tr165p7Nix6ty5s5o2bSp/f39NmDBB/v7+mjVrllq3bq3//Oc/mjNnji5evKhly5bZPD8mJkaffvqpOnfurJYtW2rmzJkqVqyYZsyYkeExOeLm5mbtDbFv3z6bdWazWZ9//rnMZrMGDRqU7rA4y0VzZm+/nhZXV1c1bNhQgwcP1i+//KKgoCB98cUXevLJJzM9p4yUN851ZuSF1/9Ojysr+xw7dqyOHDmS5f8yCqq8vb31yCOPSJJ+/vlnu/Xr1q1TXFycihYtqsaNG2dYZ061aXHmzJlMDYMEACAnERIBAPKtRx99VJK0e/fuDLfNysVZZtu19DL4888/FR0drc6dO6c50bYkFS5cWI0aNbI+vnnzpnbv3q22bdvKbDYrOTlZycnJKlasmPz9/e2CAy8vLz3++OPWxx4eHmrVqpX27t1rt6/ly5crICBArVq1SvcYnnzySUm3b3keFxcn6fZcOUOHDtXWrVvVsmVLde7cOc3nnzx5UpcuXVKlSpX04IMPpruvrCpVqpQ6deqkKVOmqEGDBpl+Xl441/fK65/V48qOfWanvn37ymQyaenSpfrtt9+syw8fPmy9o1rv3r1teuRJtyfnbtWqlcaNG5dtbR44cEBTpkyxDkdLKSgoSG+88YaSk5P1+OOPZ9jzEQCAnMJwMwBAvtWoUSOVL19eISEhunTpkkqUKJHmtpm5JXpW261cubJ1TptChQrp7bffTrfd4sWL2zyOiYnRrVu3NGHCBE2YMMFue29vb5vHjiYJLlGihKKiojI6pDS99tprWr58ufbs2aMWLVqoXLlyOnHihOLj49W8eXNNnjw53edb5szp2LFjpvY3fvx4/fXXX1muc8yYMapfv36mts0L5zqzcvv1T+1Ojutu95mdHn74YQ0aNEiTJ0/We++9p8mTJ8vLy0v//POPzGazWrRooddee83uedHR0YqMjFR0dHS2tXn9+nVNnz5d06dPV8mSJVW6dGklJSXp3Llz1smua9eubQ2aAABwBkIiAEC+ZTKZ1LlzZ02cOFFr1qxRz549HW6XmJioEydOSEr/luhZbdfFxUXly5fXwYMH9cYbb9iFQI7aTalw4cJycXHR66+/bu1RklKhQoVsHju6oL106ZLDO2NlZuJm6fbdnBYtWqQJEyZo+/btOnHihKpUqaIXX3xRnTp1ynBS5tWrV8vd3V0vvPBCuttZREVFWc9FVty4cSNT2+WVc51XX//Usnpc2bHP7NavXz9Vq1ZN8+bN08GDB3Xp0iX5+/urQ4cO6t69u1xdXXOlzWrVqumjjz7Stm3b9O+//+r48eNKSkpS0aJFFRgYqLZt26p9+/YOJ8MGACC3mAzDMJxdBAAAOSUuLk6tW7eWj4+P1q5dKxeX7BlpnZl2//nnHz377LMqVaqU1q9fn+YtyqX/uzvZunXrbJZ369ZNhQsX1syZM9OtZ9iwYVqxYoVmz55tHXKWmJio1q1bq1GjRpo0aZLN9n379tXff/+t//3vfzZD3LLT9u3b9corr+jll1/WJ598kiP7yA3Zfa6lvPP6nzlzRk888YQaN26sBQsW2K3P6nHll3MOAMD9ijmJAAD5mre3t/r166eIiAitXr06V9sdP368DMPQu+++m+HFdVo++OADbdu2Tf3799f69esVEhKi1atX6+OPP7YO67Hw8fHRp59+qp9//lkbN25Uv379dPnyZfXt29eu3T179qhWrVo5FlBI0vTp0+Xl5aX+/fvn2D5yQ06c63vl9c/qceWXcw4AwP2K/qwAgHyva9euiouLk9lszrV2586dq6CgIDVo0EDPPffcHe+jTp06Wrx4saZMmaKRI0fqxo0bKl26tBo1amQ3XKpYsWIaOXKkxo4dqxMnTqh8+fL69ttvVbVqVZvtjh8/rujoaI0YMeKO68pIXFycmjRpop49e6Y7F9S9IjvP9b3y+mf1uPLbOQcA4H7EcDMAALLJ3r17NXHiREVFRSkiIkKFCxfWsmXLrLcxz0lpDVdDznDmuc5OqYeb5ZfjAgAAd4bhZgAAZJNdu3YpJCREly9fVvPmzfXDDz9wcZ1P5ddznV+PCwAAZA49iQAAyAfoSQQAAIC7RUgEAAAAAAAAhpsBAAAAAAAgkyGRYRiKj48XnY4AAAAAAADyp0yFRDdu3NArr7yiGzdu5HQ9ACQZN2/qWo+uutajq4ybN51dDgAAAADgPsBwMwAAAAAAABASAQAAAAAAgJAIAAAAAAAAIiQCAAAAAACACIkAAAAAAAAgQiIAAAAAAACIkAgAAAAAAAAiJAIAAAAAAIAIiQAAAAAAACBCIgAAAAAAAIiQCAAAAAAAACIkAgAAAAAAgAiJAAAAAAAAIEIiAAAAAAAAiJAIAAAAAAAAIiQCAAAAAACACIkAAAAAAAAgQiIAAAAAAABIcnN2AQAAIG8wTZzq7BKcyhgy0NklAACAHGI2m/X888+rdOnSmjNnjrPLcejkyZNq27atPvroI3Xr1s0pNRASAQAAAACAe4JhGGrYsKHi4uIytb2Li4t27dqltWvX6siRIxo9enQOV3jnypcvr2effVbTpk3Tc889J29v71yvgZAIAAAAAIB8JC/3Dr7bnrvXr19Xr169bJZFRUVpyZIl8vPz03PPPWezztvbWwUKFND06dPVuHFj1alT5672n9Nef/11/fLLL5o/f77eeuutXN8/IREAAAAAALgneHt7a+BA26Bp1apVWrJkiZo3b263TpL++usvnT171imhS1b5+/urWrVqWrp0qfr27SsXl9ydSpqJqwEAAAAAwD0rPDxcklSzZk2H61esWCGTyaQnn3zSbt3UqVMVEBCgVatWOXzuTz/9pICAAM2bNy/b6s1I27ZtdfbsWW3bti3X9mlBSAQAAAAAAO5ZlpCoRo0adusMw1BISIgqVaqkIkWK2K2vX7++JGnv3r1262JiYjR58mRVqlQpVyeSrlevniRp+/btubZPC0IiAAAAAABwzzp06JDc3d1VtWpVu3XHjh3T1atXHQZI0u1AxsXFRfv377dbN3XqVEVHR+vDDz+Uu7t7ttedllq1akmSwsLCcm2fFoREAAAAAADgnnTmzBnFxMSoatWq8vDwsFt//vx5SVKJEiUcPt/b21tVqlTRoUOHlJSUZF3+zz//aNGiRWrZsqUCAwNzpvg0eHt7y9PT01p7biIkAgAAAAAA96T0hppJ0tWrVyXJ4VAzi/r16yshIUFHjhyxLhs9erRMJpOGDRuWfcVmgY+Pj6Kjo3N9v4REAAAAAADgnnTw4EFJaYdEBQoUkCQlJCSk2YZlXqJ9+/ZJktavX6/t27frlVdeUYUKFdJ8XsqeR9ktISFBBQsWzLH20+KW63sEAAAAAADIBocOHZKU9p3NfH19Jf1fjyJHLBNF79u3Tx07dtS4ceNUsmRJ9evXz2a7Hj16KCAgQMnJyVq9erUeffRRTZ48WWazWd9++62WLVum6OhoValSRUOHDlXjxo1tnnfz5k2tXr1ahQoVUp8+fdSjRw+H9ZjNZsXGxqpKlSpZeSmyBT2JAAAAAADAPSk8PFyurq6qVq2aw/VVq1aVi4uLTp48mWYbFStWlK+vr/bt26c5c+YoMjJS7777rry9ve22XbZsmXx9fbV06VINHjxYkvTtt99q/fr1+vLLL7Vy5Uo988wzeuONN3T69Gmb5xUvXlzLli1Tv379NHbsWG3dutVhPRERETKbzfL398/KS5Et6EkEAAAAAADuORcvXtTFixdVtWpV67Cy1IoUKaKAgAAdOHBAhmHIZDI53K5evXratGmT5syZo9q1a+uFF15wuF2VKlU0aNAg6+OEhAR99913WrBggWrXri1J6tWrlzZu3Kjly5dbty1fvrw1VKpUqZJCQ0M1f/58NWvWzG4flmFvlp5IuYmeRAAAAAAA4J6T0aTVFq1bt1ZsbKzD29xb1K9fX2azWQkJCRoxYkSaYVLqYW2nTp3SjRs31LNnT9WvX9/6386dO3Xq1CnrdpYAyaJu3bo6fvy4w30EBwfL1dVVLVq0SPe4cgI9iQAAAAAAwD0nsyFR586d9e2332rlypWqU6eOw23Kli0rSWrfvr11jiJHUk8mHR8fL0n67rvvVKJECZt1Xl5e1v9PHTql1avpxo0b2rBhg1q2bKnSpUunfVA5hJAIAAAAAIB8xBgy0Nkl5Ip+/frZTS7tSOnSpdW2bVutWrVK7777rk14I90ObObNmycvLy+99957WaqhcuXKcnd314ULF/Twww+nuZ1lCFnKxxUrVrTbbvXq1YqPj1evXr2yVEd2YbgZAAAAAADI19555x3Fx8dr4cKFdusWLlyo/fv36+23385y7x1vb2/17NlTo0eP1qpVq3T69Gnt379fs2bNUnBwsHW7kydPavLkyTpx4oQWLVqkdevWqXv37jZtJScna9asWWrVqpUaNWp0Zwd6l+hJBAAAAAAA8jU/Pz+NGzdO0dHRkqTz589r5cqVOnHihFauXKnGjRvrlVdeuaO2hwwZomLFimnq1Kk6e/asihYtqnr16ql169bWbTp27KgLFy6oQ4cOKliwoN5//301b97cpp3z58+rffv2eu655+78QO8SIREAAAAAAMj3nn76aev/b9myRZMmTVLRokX1zDPPaPjw4XJxSX+w1YIFCxwud3FxUe/evdW7d+80n+vu7q4xY8ZozJgxaW7z0EMPaeBA5w4VJCQCAAAAAAD3lU6dOqlTp07OLiPPYU4iAAAAAAAA0JMIAAAAAAAgp6Q1TC0voicRAAAAAAAACIkAAAAAAABASAQAAAAAAAAREgEAAAAAAECERAAAAAAAABAhEQAAAAAAAERIBAAAAAAAABESAQAAAAAAQIREAAAAAAAAECERAAAAAAAAREgEAAAAAAAAERIBAAAAAABAhEQAAAAAAAAQIREAAAAAAABESAQAAAAAAPI5s9ms9u3b64033rjjNk6ePKkaNWpo4cKF2VhZ3uLm7AIAAAAAAAAywzAMNWzYUHFxcZna3sXFRbt27dLatWt15MgRjR49+o73Xb58eT377LOaNm2annvuOXl7e99xW3kVIREAAAAAALgnXL9+Xb169bJZFhUVpSVLlsjPz0/PPfeczTpvb28VKFBA06dPV+PGjVWnTp272v/rr7+uX375RfPnz9dbb711V23lRYREAAAAAADgnuDt7a2BAwfaLFu1apWWLFmi5s2b262TpL/++ktnz57NllDH399f1apV09KlS9W3b1+5uOSvWXzy19EAAAAAAID7Snh4uCSpZs2aDtevWLFCJpNJTz75pN26qVOnKiAgQKtWrXL43J9++kkBAQGaN2+edVnbtm119uxZbdu27e6Lz2MIiQAAAAAAwD3LEhLVqFHDbp1hGAoJCVGlSpVUpEgRu/X169eXJO3du9duXUxMjCZPnqxKlSqpW7du1uX16tWTJG3fvj07ys9TGG4GAAAAAEA+YRiGlJDg7DLS5ukpk8mUrU0eOnRI7u7uqlq1qt26Y8eO6erVq3rsscccPrdevXpycXHR/v377dZNnTpV0dHRGj9+vNzd3a3La9WqJUkKCwvLpiPIOwiJAAAAAADIBwzDUPyoT3Xrn6POLiVNrv4B8hrxSbYFRWfOnFFMTIxq1KghDw8Pu/Xnz5+XJJUoUcLh8729vVWlShUdOnRISUlJ1jDon3/+0aJFi9SyZUsFBgbaPcfT09Padn7CcDMAAAAAAPKLbO6lk9elN9RMkq5evSpJDoeaWdSvX18JCQk6cuSIddno0aNlMpk0bNgwh8/x8fFRdHT0HVadd9GTCAAAAACAfMBkMslrxCf31XCzgwcPSko7JCpQoIAkKSGd16R+/fpavHix9u3bp1q1amn9+vXavn27evfurQoVKjh8TkJCggoWLHh3xedBhEQAAACSTBOnOrsEOIkxxP52yQBwrzKZTNL/D0buB4cOHZKU9p3NfH19Jf1fjyJHLBNR79u3Tx07dtS4ceNUsmRJ9evXz+H2ZrNZsbGxqlKlyp0Xnkcx3AwAAAAAANyTwsPD5erqqmrVqjlcX7VqVbm4uOjkyZNptlGxYkX5+vpq3759mjNnjiIjI/Xuu+/K29vb4fYREREym83y9/fPlmPISwiJAAAAAADAPefixYu6ePGiKlWqZB1WllqRIkUUEBCgAwcO3L7zWxrq1aunEydOaM6cOapdu7ZeeOGFNLfdt2+fJKlx48Z3dwB5ECERAAAAAAC452Q0abVF69atFRsb6/A29xb169eX2WxWQkKCRowYke68ScHBwXJ1dVWLFi3uqO68jJAIAAAAAADcczIbEnXu3Fmurq5auXJlmtuULVtWktS+fXvrHEWO3LhxQxs2bFDLli1VunTprBedxxESAQAAAACAe06/fv105MgR9erVK93tSpcurbZt22rVqlWKj4+3W28YhubNmycvLy+999576ba1evVqxcfHZ7jPexUhEQAAAAAAyNfeeecdxcfHa+HChXbrFi5cqP379+vtt99Ot3dQcnKyZs2apVatWqlRo0Y5Wa7TuDm7AAAAAAAAgJzk5+encePGKTo6WpJ0/vx5rVy5UidOnNDKlSvVuHFjvfLKK+m2cf78ebVv317PPfdcbpTsFIREAAAAAAAg33v66aet/79lyxZNmjRJRYsW1TPPPKPhw4fLxSX9wVYPPfSQBg4cmNNlOhUhEQAAAAAAuK906tRJnTp1cnYZeQ5zEgEAAAAAAICQCAAAAAAAAIREAAAAAAAAECERAAAAAAAAREgEAAAAAAAAERIBAAAAAABAhEQAAAAAAAAQIREAAAAAAABESAQAAAAAAAAREgEAAAAAAECERAAAAAAAABAhEQAAAAAAAERIBAAAAAAAABESAQAAAAAAQIREAAAAAAAAECERAAAAAAAAREgEAAAAAAAAERIBAAAAAABAhEQAAAAAAAAQIREAAAAAAABESAQAAAAAAAAREgEAAAAAAECERAAAAAAAABAhEQAAAAAAAERIBAAAAAAAABESAQAAAAAAQIREAAAAAAAAECERAAAAAAAAREgEAAAAAAAAERIBAAAAAABAhEQAAAAAAAAQIREAAAAAAABESAQAAAAAAAAREgEAAAAAAECERAAAAAAAABAhEQAAAAAAAERIBAAAAAAAABESAQAAAAAAQIREAAAAAAAAECERAAAAAAAAREgEAAAAAAAAERIBAAAAAABAhEQAAAAAAAAQIREAAAAAAABESAQAAAAAAAAREgEAAAAAAECERAAAAAAAABAhEQAAAAAAAERIBAAAAAAAABESAQAAAAAAQIREAAAAAAAAECERAAAAAAAAREgEAAAAAAAAERIBAAAAAABAhEQAAAAAAACQ5ObsAgAAyEtME6c6uwQAAADAKehJBAAAAAAAAEIiAAAAAAAAEBIBAAAAAABAhEQAAAAAAAAQIREAAAAAAABESAQAAAAAAAAREgEAAAAAAECERAAAAAAAABAhEQAAAAAAAERIBAAAAAAAABESAQAAAAAAQIREAAAAAAAAECERAAAAAAAAREgEAAAAAAAAERIBAAAAAABAkpuzCwAA5D2miVOdXQIAAACAXEZPIgAAAAAAABASAQAAAAAAgJAIAAAAAAAAIiQCAAAAAACACIkAAAAAAAAgQiIAAAAAAACIkAgAAAAAAAAiJAIAAAAAAIAIiQAAAAAAACBCIgAAAAAAAIiQCAAAAAAAACIkAgAAAAAAgAiJAAAAAAAAIEIiAAAAAAAAiJAIAAAAAAAAIiQCAAAAAACACIkAAAAAAAAgQiIAAAAAAACIkAgAAAAAAAAiJAIAAAAAAIAIiQAAAAAAACBCIgAAAAAAAIiQCAAAAAAAACIkAgAAAAAAgAiJAAAAAAAAIEIiAAAAAAAAiJAIAAAAAAAAIiQCAAAAAACACIkAAAAAAAAgQiIAAAAAAACIkAgAAAAAAAAiJAIAAAAAAIAIiQAAAAAAACBCIgAAAAAAAIiQCAAAAAAAACIkAgAAAAAAgCQ3ZxcAAAAAOJNp4lRnl+A0xpCBzi4BAJCH0JMIAAAAAAAAhEQAAAAAAAAgJAIAAAAAAIAIiQAAAAAAACBCIgAAAAAAAIiQCAAAAAAAACIkAgAAAAAAgAiJAAAAAAAAIEIiAAAAAAAAiJAIAAAAAAAAIiQCAAAAAACACIkAAAAAAAAgQiIAAAAAAACIkAgAAAAAAAAiJAIAAAAAAIAIiQAAAAAAACBCIgAAAAAAAIiQCAAAAAAAACIkAgAAAAAAgAiJAAAAAAAAIEIiAAAAAAAAiJAIAAAAAAAAIiQCAAAAAACACIkAAAAAAAAgQiIAAAAAAACIkAgAAAAAAAAiJAIAAAAAAIAIiQAAAAAAACBCIgAAAAAAAIiQCAAAAAAAACIkAgAAAAAAgAiJAAAAAAAAIEIiAAAAAAAAiJAIAAAAAAAAIiQCAAAAAACACIkAAAAAAAAgQiIAAAAAAACIkAgAAAAAAAAiJAIAAAAAAIAIiQAAAAAAACBCIgAAAAAAAIiQCAAAAAAAACIkAgAAAAAAgAiJAAAAAAAAIEIiAAAAAAAAiJAIAAAAAAAAIiQCAAAAAACACIkAAAAAAAAgQiIAAAAAAACIkAgAAAAAAAAiJAIAAAAAAIAIiQAAAAAAACBCIgAAAAAAAIiQCAAAAAAAACIkAgAAAAAAgAiJAAAAAAAAIEIiAAAAAAAAiJAIAAAAAAAAIiQCAAAAAACACIkAAAAAAAAgQiIAAAAAAACIkAgAAAAAAAAiJAIAAAAAAIAIiQAAAAAAACBCIgAAAAAAAIiQCAAAAAAAACIkAgAAAAAAgAiJAAAAAAAAIEIiAAAAAAAAiJAIAAAAAAAAIiQCAAAAAACACIkAAAAAAAAgQiIAAAAAAACIkAgAAAAAAAAiJAIAAAAAAIAIiQAAAAAAACBCIgAAAAAAAIiQCAAAAAAAACIkAgAAAAAAgAiJAAAAAAAAIMnN2QUAAAAAcA7TxKnOLsGpjCEDnV0CAOQp9CQCAAAAAAAAIREAAAAAAAAIiQAAAAAAACBCIgAAAAAAAIiQCAAAAAAAACIkAgAAAAAAgAiJAAAAAAAAIEIiAAAAAAAAiJAIAAAAAAAAIiQCAAAAAACACIkAAAAAAAAgQiIAAAAAAACIkAgAAAAAAAAiJAIAAAAAAIAIiQAAAAAAACBCIgAAAAAAAIiQCAAAAAAAACIkAgAAAAAAgAiJAAAAAAAAIEIiAAAAAAAAiJAIAAAAAAAAIiQCAAAAAACACIkAAAAAAAAgQiIAAAAAAACIkAgAAAAAAAAiJAIAAAAAAIAkN2cXACDvMk2c6uwSnMYYMtDZJQAAAABArqInEQAAAAAAAAiJAAAAAAAAQEgEAAAAAAAAMScRADh0P8/HBAAAAOD+RE8iAAAAAAAAEBIBAAAAAACAkAgAAAAAAAAiJAIAAAAAAIAIiQAAAAAAACBCIgAAAAAAAIiQCAAAAAAAACIkAgAAAAAAgAiJAAAAAAAAIEIiAAAAAAAAiJAIAAAAAAAAIiQCAAAAAACACIkAAAAAAAAgQiIAAAAAAACIkAgAAAAAAAAiJAIAAAAAAIAIiQAAAAAAACBCIgAAAAAAAIiQCAAAAAAAACIkAgAAAAAAgAiJAAAAAAAAIEIiAAAAAAAAiJAIAAAAAAAAIiQCAAAAAACACIkAAAAAAAAgQiIAAAAAAACIkAgAAAAAAAAiJAIAAAAAAIAIiQAAAAAAACBCIgAAAAAAAIiQCAAAAAAAACIkAgAAAAAAgAiJAAAAAAAAIEIiAAAAAAAAiJAIAAAAAAAAIiQCAAAAAACACIkAAAAAAAAgQiIAAAAAAACIkAgAAAAAAACS3JxdAAAAAAA4g2niVGeX4DTGkIHOLgFAHkRPIgAAAAAAABASAQAAAAAAgJAIAAAAAAAAIiQCAAAAAACACIkAAAAAAAAgQiIAAAAAAACIkAgAAAAAAAAiJAIAAAAAAIAIiQAAAAAAACBCIgAAAAAAAIiQCAAAAAAAACIkAgAAAAAAgAiJAAAAAAAAIEIiAAAAAAAAiJAIAAAAAAAAIiQCAAAAAACACIkAAAAAAAAgQiIAAAAAAACIkAgAAAAAAAAiJAIAAAAAAIAIiQAAAAAAACBCIgAAAAAAAIiQCAAAAAAAACIkAgAAAAAAgAiJAAAAAAAAIEIiAAAAAAAAiJAIAAAAAAAAIiQCAAAAAACACIkAAAAAAAAgQiIAAAAAAACIkAgAAAAAAAAiJAIAAAAAAIAIiQAAAAAAACBCIgAAAAAAAIiQCAAAAAAAACIkAgAAAAAAgAiJAAAAAAAAIEIiAAAAAAAAiJAIAAAAAAAAktycXQAAAAAAIHeZJk51dglOYwwZ6OwSgDyLnkQAAAAAAAAgJAIAAAAAAAAhEQAAAAAAAERIBAAAAAAAABESAQAAAAAAQNzdDJnAnQ8AAAAAAMj/6EkEAAAAAAAAQiIAAAAAAAAw3AxIl7OG2nnduqVz////vb+ZqXhXV6fUAQAAAAC4f9CTCAAAAAAAAIREAAAAAAAAICQCAAAAAACACIkAAAAAAAAgQiIAAAAAAACIkAgAAAAAAAAiJAIAAAAAAIAIiQAAAAAAACBCIgAAAAAAAIiQCAAAAAAAACIkAgAAAAAAgAiJAAAAAAAAIEIiAAAAAAAAiJAIAAAAAAAAktycXQAAAAAAALnFNHGqs0twKmPIQGeXgDyMnkQAAAAAAADIWk8inymzlOx2/3U+ImkFAAAAAAD5HT2JAAAAAAAAwJxEmXG/j1kFAAAAAAD5Hz2JAAAAAAAAQEgEAAAAAACATA43Mwzj9sa3knO0GAC3uZlv6UaKz53b//9/AAAAALgb7mO/dnYJThPzdh9nl+B0BQsWlMlkSnO9yTAyvvq8fPmy+vbtm62FAQAAAAAAIPf88MMP8vLySnN9pkIis9ms6OhoFShQIN3EKadcv35dgYGBCgoKUqFChXJ9/3Aezv39ifN+/+Lc37849/cvzv39i3N//+Lc3784986XUU+iTA03c3FxUfHixbOtqKwym80ym80qWLBguokX8h/O/f2J837/4tzfvzj39y/O/f2Lc3//4tzfvzj3eR8TVwMAAAAAAICQCAAAAAAAAPdISOTh4aEBAwbIw8PD2aUgl3Hu70+c9/sX5/7+xbm/f3Hu71+c+/sX5/7+xbnP+zI1cTUAAAAAAADyt3uiJxEAAAAAAAByFiERAAAAAAAACIkAAAAAAABASAQAAAAAAAAREgEAAAAAAECSm7MLSM++ffs0depU7dmzR0lJSapSpYpeeeUVPfvss84uDTnkwoULWrt2rYKCgnT8+HFdunRJPj4+atCggXr37q26des6u0Tkkjlz5mjixImSpMWLF6tevXrOLQg57o8//tCPP/6o8PBw3bhxQyVKlFC9evX0/vvv68EHH3R2ecgBhmHojz/+0IIFC3TixAnFxsbqgQceUJMmTfTGG2/Iz8/P2SXiLv3666/avXu3Dhw4oKNHjyopKUljxoxRhw4dHG4fFxenqVOn6vfff9fFixdVsmRJtWnTRgMHDpS3t3cuV4+7kdlzn5SUpL/++kt///239u7dq3PnzslkMqlKlSp6/vnn9dJLL8nV1dVJR4GsyupnPqXTp0+rffv2io+PV5cuXfT555/nQsXILndy7k+fPq1Zs2YpODhYFy9eVJEiRVS5cmW9/PLLatu2bS5Wj5TybEi0Y8cOvf7663J3d1e7du1UuHBh/f777xoyZIgiIyPVt29fZ5eIHLBgwQLNmTNH5cqVU7NmzVS8eHGdPHlSGzZs0IYNGzRp0iQ9/fTTzi4TOezYsWOaMmWKvLy8FB8f7+xykMMMw9Ann3yixYsXq1y5cnr66adVqFAhRUVFaefOnYqMjCQkyqfGjRun77//XiVLltQTTzwhb29vHT58WEuWLNFvv/2mn376Sf7+/s4uE3fhm2++UWRkpHx9fVWqVClFRkamuW18fLy6d++uQ4cO6dFHH1W7du10+PBhzZs3Tzt27NCPP/4oLy+vXKwedyOz5/7UqVN6++23VahQITVt2lStWrVSbGys/v77b33++efavHmzZsyYIZPJlMtHgDuRlc98SoZhaPjw4TlcHXJSVs99cHCw+vfvL0lq2bKl/Pz8FBMToyNHjmjbtm2ERM5k5EFJSUlG69atjVq1ahkHDx60Lo+NjTXatWtn1KhRwzhx4oTzCkSOWb9+vbFz50675Tt37jRq1qxpNG7c2EhISHBCZcgtycnJRseOHY1OnToZQ4YMMfz9/Y2wsDBnl4Uc9MMPPxj+/v7GZ599ZiQnJ9utT0pKckJVyGlRUVFGtWrVjJYtWxqxsbE2677//nvD39/fGDZsmJOqQ3YJDg42zpw5YxiGYcyaNcvw9/c3li1b5nDbb775xvD39zfGjx/vcPk333yT4/Ui+2T23J8/f95YuHChER8fb7P8+vXrRocOHQx/f39jzZo1uVIz7l5WPvMp/fDDD0aNGjWs3/8jR47M6VKRzbJy7s+ePWs0aNDAaNOmjREZGWm3nr/9nCtPzkm0fft2nTp1Ss8884xq1KhhXe7t7a233npLycnJWr58uRMrRE5p06aNGjZsaLe8YcOGatKkia5evaojR444oTLkljlz5ujw4cP68ssv6V5+H7h586amT58uPz8/DR8+3OE5d3PLs51ecRciIyNlNpvVoEEDu2FELVq0kCRduXLFCZUhOzVr1kxly5bNcDvDMLR06VJ5eXlZf1m26NOnj3x8fPTzzz/LMIycKhXZLLPnvnTp0nr55ZdVsGBBm+VeXl569dVXJUk7d+7MkRqR/TJ73lM6efKkvvrqK/Xu3VvVq1fPocqQ07Jy7mfOnKm4uDh9+umnKlOmjN16/vZzrjz56oeEhEiSmjdvbrfu0UcftdkG9w/LlwVfGvnX0aNHNW3aNPXr109Vq1Z1djnIBcHBwbp69apeeOEFmc1m/f7774qIiFDhwoXVrFkzlS9f3tklIoeUL19e7u7uCg0NVVxcnE1QtGnTJklS06ZNnVUecllERISioqLUvHlzuyFlnp6eatiwof7880+dPHlSFSpUcE6RyHWWv/n40Sj/MpvN+vDDD1WmTBn1799fYWFhzi4JOcwwDK1bt05FixbVI488ogMHDmjnzp0ym82qXr26mjZtKheXPNmX5b6RJ6+2IyIiJMnhxYGPj498fX118uTJXK4KznT27Flt3bpVJUuWZH6KfCo5OVnDhg1T5cqV9eabbzq7HOSSAwcOSLp9AdC+fXudOHHCus7FxUW9evXSBx984KzykIN8fX01ePBgjR8/Xk8//bRatWqlQoUK6ejRo9q2bZu6dOmi7t27O7tM5BLL33VpBUCWvwkJie4vy5Ytk+T4h2PkDz/88IPCwsL0448/ysPDw9nlIBecOXNGV69eVe3atfXJJ5/op59+sllfo0YNzZgxQw888ICTKkSeDIni4uIkSYULF3a43tvbW+fPn8/NkuBESUlJGjp0qBITEzVkyBB+TcqnZs6cqSNHjmjJkiVyd3d3djnIJZcvX5Ykff/996pRo4aWLl2qypUr69ChQxo5cqTmzp0rPz8/vfzyy06uFDnh9ddfV6lSpfTxxx9r0aJF1uX169dX+/bt+S64j8TGxkpSmncwsyy3bIf8b/HixQoKClLTpk31+OOPO7sc5IATJ05o8uTJ6tmzp+rXr+/scpBLLH/7hYeH69ixYxozZoyeeOIJxcbGatasWVqyZInefvttLVmyxMmV3r/ox4U8zWw2a/jw4dq5c6defPFFPf/8884uCTng8OHDmjlzpl577TXVrFnT2eUgF1nmF3F3d9f06dNVp04dFSpUSA0bNtSUKVPk4uKi77//3slVIqd8++23+vDDD9WnTx9t2rTJ+mvyrVu31LNnT/3+++/OLhGAE2zcuFGjRo1S2bJlNWHCBGeXgxxgGWZWqlQpvfPOO84uB7nIbDZLkm7duqVBgwapQ4cO8vHx0UMPPaRRo0apbt262rt3r3bt2uXkSu9feTIkyujXori4uDR7GSH/MAxDI0aM0MqVK9W+fXt99tlnzi4JOeSDDz6Qn5+fBg4c6OxSkMss3/e1atVS6dKlbdZVrVpVfn5+OnXqlK5du+aM8pCDtm3bpm+++UbdunVT37599cADD8jLy0sPP/ywZs2aJU9PT40ZM8bZZSKXWP6us/QmTy2jXubIPzZv3qyBAweqePHi+uGHH1SqVClnl4QcMH/+fO3Zs0ejR4+2m7Qc+VvK7/EnnnjCbn3Lli0l/d+UBMh9eXK4mWWs+cmTJ1WrVi2bdTExMYqOjqZLYj5nNpv10Ucfafny5XrmmWc0duxYJjDLxw4fPixJql27tsP1Xbp0kSRNnz5drVu3zrW6kPMqVaokKe0LP8vymzdvqkiRIrlWF3KeZXLqJk2a2K0rVqyYAgICFBYWpitXrqhYsWK5XR5ymWXOIcu8lKlZ5ixiMvv8LSgoSAMGDJCvr6/mz58vPz8/Z5eEHHL48GEZhqGePXs6XL948WItXrxYTzzxhL799ttcrg45qXz58nJ1ddWtW7cc/v1n+XsvISEht0vD/5cnQ6JGjRpp1qxZ2rJli9q1a2ezLjg4WJLUuHFjZ5SGXJAyIHr66ac1fvx45iHK5zp16uRw+a5duxQREaFWrVqpWLFiWb6lKvI+S0Bw/Phxu3VJSUk6deqUvLy8CAnyoaSkJElp3+bespyJTO8PFSpUUKlSpRQaGqr4+HibO5wlJCRo165dKlWqFCFRPhYUFKT+/fvLx8dH8+fP51znc40aNXL49/3Fixe1adMmVapUSQ0aNFCNGjWcUB1ykoeHh+rXr69du3bp33//VcOGDW3W//vvv5LE3/1OlCdDokceeUR+fn767bff1LNnT1WvXl3S7a7G3377rdzc3PTCCy84uUrkhJQB0VNPPaUJEyYQEN0HvvjiC4fLhw0bpoiICPXp00f16tXL3aKQK8qVK6fmzZtry5YtWrp0qTp37mxdN3v2bF27dk3t27e33gYZ+UeDBg30v//9T/PmzdOTTz5p82viihUrdPLkSdWsWTPNiYyRv5hMJnXu3FnTp0/X9OnT9f7771vXzZo1SzExMerfv79MJpMTq0ROSR0QcQe7/K9jx47q2LGj3fIdO3Zo06ZNatSokT7//HMnVIbc0LVrV+3atUvTpk3T7NmzrT8IHTt2TCtWrFChQoX02GOPObnK+1ee/Kvbzc1No0ePVu/evfXyyy/rmWeekbe3t37//XedOXNG77zzjipWrOjsMpEDpk+fruXLl8vLy0sVKlTQjBkz7LZp3bq1NTgEcO/75JNP9NJLL2nEiBHasGGDKlWqpPDwcG3fvl1ly5bV0KFDnV0icsBTTz2ln376SSEhIWrTpo1atWqlIkWK6MiRIwoODpaHh4eGDx/u7DJxl5YuXardu3dLko4ePWpdFhISIun2v+mWYcS9e/fWX3/9pe+++06HDh1SzZo1dfjwYQUFBal69erq3bu3cw4CdySz5/7YsWPq37+/EhMT1bhxY61evdqurbJly6pDhw65VzzuWFY+88hfsnLu27Vrp99//13r169X+/bt1bx5c8XFxWn9+vVKSEjQuHHj5OPj45wDQd4MiSSpadOm+vHHHzVlyhStXbtWSUlJqlKligYNGqT27ds7uzzkkMjISElSfHy8Zs6c6XCbsmXLEhIB+Ui5cuW0bNkyTZkyRZs3b1ZwcLBKlCihbt26qX///ipevLizS0QOcHV11X//+1/98MMPWrt2rVavXq2kpCQVL15czzzzjPr06SN/f39nl4m7tHv3bq1YscJmWWhoqEJDQyXd/jfdctHg5eWlBQsWaNq0aVq/fr1CQkJUokQJ9erVSwMGDLAZgoa8L7Pn/tKlS0pMTJQkhwGRdHuaCUKie0NWPvPIX7Jy7k0mk7766iv973//088//6zFixfLw8NDDRo0UJ8+fZhaxslMhuX+wwAAAAAAALhvcbsoAAAAAAAAEBIBAAAAAACAkAgAAAAAAAAiJAIAAAAAAIAIiQAAAAAAACBCIgAAAAAAAIiQCAAAAAAAACIkAgAAAAAAgAiJAAAAAAAAIEIiAAAAAAAAiJAIAAAAAAAAIiQCAAAAAACApP8Hl6skdO8c6t0AAAAASUVORK5CYII=",
      "text/plain": [
       "<Figure size 1150x500 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "# plot\n",
    "plt.hist(\n",
    "    Tyr,\n",
    "    np.arange(19),\n",
    "    align='left',\n",
    "    label=r'$T(y_\\mathrm{rep})$',\n",
    "    color='C0',\n",
    ")\n",
    "plt.axvline(Ty, color='C2', label='$T(y)$')\n",
    "plt.xlim((-0.5, 17.5))\n",
    "plt.yticks([])\n",
    "plt.title(\n",
    "    'Binomial example - number of changes? \\n'\n",
    "    r'$\\operatorname{Pr}(T(y_\\mathrm{rep},\\theta) \\leq T(y,\\theta)|y) = 0.03$',\n",
    ")\n",
    "plt.legend();"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "bda_py_demos",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.12.11"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 0
}
